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Abstract 

We study the zero-temperature phase diagram of a one-dimensional array of 
QED cavities where, besides the single-photon hopping, an additional coupling 
between neighboring cavities is mediated by an AMype four-level system. By 
varying the relative strength of the various couplings, the array is shown to exhibit 
a variety of quantum phases including a polaritonic Mott insulator, a density-wave 
and a superfluid phase. Our results have been obtained by means of numerical 
density-matrix renormalization group calculations. The phase diagram was 
obtained by analyzing the energy gaps for the polaritons, as well as through a 
study of two-point correlation functions. 

Keywords: Cavity quantum electrodynamics; Strongly correlated polariton 
systems; Quantum phase transitions 

PACS Codes: 42.50.Pq; 85.25.Cp; 64.70.Tg 


1 Introduction 

The recent impressive advances in the field of quantum simulators allowed to probe the 
many-body physics of strongly correlated systems at the level of the single quantum ob¬ 
ject. At present cold atoms trapped in optical lattices can be considered among the most 
promising examples of quantum simulators. By means of ultracold atomic and molecular 
gases, it is nowadays possible to reach a degree of control and accuracy in engineering the 
dynamics of many-body systems that were unimaginable in previously. As a consequence, 
the coherent quantum dynamics emerging from carefully tailored microscopic Hamiltoni¬ 
ans can now be tested experimentally [1], It has been possible, just to recall one example, 
to implement the Bose-Hubbard (BH) model [2,3] and to detect its zero-temperature su¬ 
perfluid (SF) to Mott insulator (MI) quantum phase transition [4]. Other models involving 
spinor gases, Fermi systems, Bose-Fermi mixtures, or dipolar gases have been also devised 
and realized, providing an even richer phase diagram (see for example the review [5]). We 
mention the stabilization of density-wave (DW) phases for bosons, as well as more pe¬ 
culiar topological or supersolid orderings, which can arise in the presence of finite-range 
interactions [6]. 

More recently a novel kind of many-body quantum simulator has been introduced, based 
on the idea to use single photons as quantum objects. Since photons hardly interact in open 
space, the most natural way to significantly increase their interactions is to trap them into an 
optical QED cavity, and couple the field with atoms/molecules inside it in order to create an 
optical nonlinearity. If the nonlinearity is sufficiently large, the so called photon blockade 
sets in [7, 8], namely, the presence of a single photon inside a cavity prevents a second one 
to enter it. In the rotating-wave approximation, the simplest light-matter interaction scheme 
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of this type can be accurately described by the Jaynes-Cummings (JC) model. By arranging 
an array of cavities coupled through the photon hopping, such to generate a competition 
between the hopping and the on-site nonlinearities, one can devise a setup that is well 
described by the so called Jaynes-Cummings-Hubbard (JCH) model [9-11]. 

In many respects, if one ignores dissipation, the physics emerging from the JCH Hamil¬ 
tonian resembles, at low-energies, that of an effective BH model. Probably the main differ¬ 
ence between the two systems is that, instead of having neutral bosons as building blocks 
of the model, in the JCH Hamiltonian one has to think in terms of polaritons, i.e ., combined 
photonic/atomic excitations. Many different works already addressed the JCH equilibrium 
phase diagram with analytical, as well as numerical methods, leading to a fairly complete 
theoretical understanding of the nature and the location of the emerging phases and phase 
transitions in terms of the parameters governing the system (the field has been recently 
reviewed in, e.g.. Refs. [12-15]). 

Additional interest in cavity arrays comes from the fact that these systems can be nat¬ 
urally considered as open-system quantum simulators. Some related features have been 
recently explored [16-24]. In the following we will not touch on this and consider only the 
“equilibrium” phase diagram. 

This intense activity has been very recently boosted by the first experiments on QED 
cavity arrays [25-27]. As of today, the most concrete possibility to realize controllable and 
scalable quantum simulators with cavity arrays involves circuit-QED cavities [28-30]. 

So far the coupling between cavities has been mostly considered through photon hop¬ 
ping. Only few works started addressing more general schemes, where the cavity coupling 
can be induced also by means of non-linear elements [23, 24,31, 32]. Such configurations 
include cross-Kerr interactions and/or correlated hopping terms, which lead to general¬ 
izations of the JCH model in a way similar to the extended BH (EBH) Hamiltonian for 
atoms with large dipole momentum loaded in optical lattices [33]. The underlying phys¬ 
ical model is believed to possess a much richer structure, with the emergence of exotic 
phases of correlated polaritons. It is particularly interesting to address these schemes in 
one-dimensional (ID) systems, where interactions become crucial to stabilize exotic phases 
of matter [33-37], These notably include a series of nontrivial density-wave (DW) states, 
which can arise in the strong coupling regime [38], as well as supersolidity and phase- 
separation effects [39,40]. Extension to consider also counter-rotating terms in the ultra¬ 
strong coupling regime, thus leading to the so called Rabi-Hubbard model [41], have been 
investigated [42]. However we are not aware of numerical investigations of coupled cavity 
models beyond the JCH and Rabi-Hubbard model. 

In all such situations, non-perturbative, either numerical or analytical calculations are 
necessary. Here the density-matrix renormalization group (DMRG) algorithm [43,44] has 
been employed to work out the quantitative zero-temperature phase diagram of the JCH 
model [45-47]. This is a particularly efficient method for the statics of ID many-body 
problems. Its key strategy consists in constructing a portion of the system (called block) 
and then recursively enlarge it. At each step, the basis of the corresponding Hamiltonian is 
truncated to a given value m, so that one can manage the Hamiltonian in an effective Hilbert 
space of fixed dimensions, as the physical system grows. This truncation is performed by 
retaining the eigenstates corresponding to the m highest eigenvalues of the reduced density 
matrix of the block. 

The aim of this paper is to quantitatively study a generalization of the JCH Hamiltonian, 
aimed at taking into account an effective nearest-neighbor nonlinearity between cavities 
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mediated by an IV- type four-level system as discussed for two cavities in Ref. [48]. The 
presence of this coupling leads to an effective cross-Kerr non-linearity. An analysis at the 
mean-field level of a dissipative open EBH as an effective model for nonlinearly coupled 
cavities has been performed, unveiling the emergence of novel photon crystal and super¬ 
solid phases [23,24], Here we do not resort to the effective EBH model and analyze the full 
model as introduced in [48], Using the DMRG algorithm, we work out the ID ground-state 
phase diagram. We show that a physics similar to the EBH model appears, with a rich phase 
diagram including gapless SF, as well as MI and DW phases of polaritons. We postpone 
the analysis of the interplay of driving and dissipation to a future work. 

The paper is organized as follows. In the next two sections we introduce the model of 
coupled cavities of our interest (Sec. 2) and the quantities we are going to address, namely 
the energy gaps, and the staggered number-number correlations (Sec. 3). In Sec. 4 we 
discuss the zero-temperature equilibrium phase diagram, focusing on the MI/SF boundary 
and on the boundary separating the DW from the other phases. Finally, in Sec. 5 we draw 
our conclusions. 

2 The model 

Let us consider a ID array of QED cavities, where photons can hop between neighboring 
cavities. Moreover two adjacent resonators are also nonlinearly coupled to each other via a 
A'-type four-level system, as shown in Fig. 1(a). For the sake of clarity in our description, 
we shall divide the ID array into coupled effective sites composed of a cavity and an atom. 
The four levels are denoted by {U) }/= i ...4, and are depicted in Fig. 1(b). An external laser 



(b) |4) 




Figure 1 Scheme for nonlinearly coupled QED cavities, (a) An array of QED cavities nonlinearly 
coupled by A-type atoms. The photon hopping between nearest-neighbor cavities has a strength t. 

Each effective site is composed of a cavity and an atom (dashed box), (b) Level structure of the AMype 
atoms. The transition |1) <-* |3> is resonantly coupled to the cavity mode of its own site with strength gi, 
while |2) <-> |4> is coupled to the cavity mode of its right nearest-neighbor site with strength g 2 , and has a 
detuning A. The transition |2) <-» |3> is in resonance with an external laser field of strength £2. 
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with frequency Li resonantly drives the transition |3) |2). The transition |1) <-> |3) is 

resonantly coupled to the cavity mode of the same site with strength g i, while the transition 
|2) |4) couples to the cavity mode of its right nearest-neighbor site with strength gi, and 

a detuning A. 

The use of such /V-type atom for generating large Kerr nonlinearity has been extensively 
studied in the literature [7, 8,49], however the vast majority of the scenarios only focused 
on a single-mode cavity. Our work is inspired by the idea of Ref. [48], where the cross-Kerr 
nonlinearity is generated between two different and neighboring cavities, in circuit-QED 
systems. In practice, we use the unbalanced couplings of atomic transition |1) |3) with 

left cavity mode, and |2) «-» |4) with right cavity mode respectively, in order to generate 
the local (gi) and nonlocal (g 2 ) nonlinearities of our many-body system. This kind of four- 
level artificial molecule can be realized using two Josephson transmon qubits coupled by a 
superconducting quantum interference device. 

Using the interaction picture and in the rotating-wave approximation, the system Hamil¬ 
tonian reads 
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where cr" m = | m) (n \, ( m , n = 1,2,3,4), and a (af is the annihilation (creation) operator of 
the cavity mode. The subscripts denote the site position along the ID chain. The first three 
terms in the r.h.s. of Eq. (1) describe the local terms and the nonlinearities on each site. 
Inside the latter brackets, the first term is the photon hopping, while the second term de¬ 
scribes the coupling of the atom to its right neighboring cavity, which generates an effective 
nonlocal cross-Kerr nonlinearity between the two cavities. 

Hereafter we concentrate on the ID model in Eq. (1) at zero temperature, specifically 
addressing the case without dissipation with DMRG. Let us also fix the Hamiltonian quan¬ 
tities in units of D, set h = 1, and work with open boundary conditions. We recall that, 
in the presence of dissipation, the problem becomes much more difficult to be handled 
numerically 

For the system we are considering here, in the strong coupling regime atoms and photons 
cannot be considered as two separate entities. It is thus natural to investigate the phase 
diagram in terms of combined atomic/photonic modes, named polaritons. The polaritonic 
number operator on each site i, representing the number of local excitations, is defined as 
nf = 2crf + cr 33 + of 2 + a\ a ;. For the closed system described by the Hamiltonian (1), the 
total number /V po1 = nf of such polaritons is a conserved quantity. In the following we 
work in the canonical ensemble for polaritons, and focus on the integer filling situation. 


™It is however possible to address the effect of dissipation with a DMRG approach in 
a ID chain, when this is described by a master equation within the Lindblad formal¬ 
ism. In the language of tensor networks, one has to generalize the matrix-product-state 
ansatz to a matrix-product-density-operator ansatz for mixed states, as originally proposed 
in Refs. [50,51]. The computational complexity is greater than for static computations, and 
is eventually related to the amount of entanglement in the steady state. 
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3 Energy gaps and correlation functions 

The different nature of the various phases is sensitive to a number of properties which we 
are going to focus on. Here we are going to study quantities that resemble those character¬ 
izing the various phases of the EBH model [36]. 

First of all, the ground-state energy gap is an important indicator which characterizes the 
presence or absence of criticality in the model. In particular, in the critical SF phase, the 
charge gap vanishes in the thermodynamic limit. On the other side, in the insulating MI and 
DW phases, such gap remains finite. In order to make connection with a similar notation 
in the EBH model, below we introduce the so called charge and neutral gaps referring 
respectively to the gaps corresponding to adding one extra particle (“charge” sector) or 
remaining with the same number of particles (“neutral” sector). We stress however that in 
the present model the excitation carry no real charge. This has to be understood only as a 
convention. 

The charge gap is defined as 


A E c = A E + - A E~ , (2) 

where, in the canonical ensemble, A E + (A/: ) denotes the extra energy needed to add (re¬ 
move) one particle, i.e. one polariton, in the system. In the specific, focusing on the unit 
filling, A E + — E / + \ - E r and A/: = E L - E r where E t is the ground-state energy 

per site of an L-sites cavity-array with exactly L excitations, and E lA \ (E r \) is the corre¬ 
sponding energy per site with one excitation more (less). It is therefore possible to extrap¬ 
olate A E c by running three different DMRG simulations with fixed number of polaritons 
N poi =L- 1, L, L+ 1 [52,53], 

While the charge gap is able to detect particle-hole excitations, in some circumstances it 
is possible that the dominant low-energy excitations are of a different type. Their presence 
can be detected only by the so called neutral gap at a fixed number of particles, 

A E„ = E) - E l , (3) 

where, again working in the canonical ensemble, E l L denotes the first excited energy per 
site of an L-site system with L excitations. 

In the following, we also focus on the analysis of the staggered diagonal order for the 
polaritons, in order to distinguish the DW from the other phases. We do this by investigating 
the two-point correlation function 

CD W (r) = (-1) r <bnf%^>, (4) 

where 5n p ° l = n? 01 —h denotes the polariton fluctuation from the average filling h. The order 
parameter identifying the DW phase is thus given by: (?dw = hm MOO CdwM- A finite 
value of <9 dw indicates a tendency to establish, in the thermodynamic limit, a staggered 
occupation of the polaritons. On the other side, in the MI as well as the SF phases, CdwM 
vanishes exponentially with increasing distance r. 

4 Phase diagram 

The zero-temperature phase diagram of model (1), at unit polariton filling h = 1 and in 
the g 2 - t plane, is summarized in Fig. 2. We observe that three different phases can be 
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Figure 2 Phase diagram. Ground-state phase diagram for a 1 D system of coupled cavities described by 
Hamiltonian (1) in the g 2 - t parameter space. Here and in the subsequent figures we choose 
gi/fi = 0.8, A/£l = -2. The various symbols correspond to points belonging to different phases 
(triangles = Ml, circles = DW, squares = SF). Filled points lie close to a phase transition, and are used to 
draw the interpolating curves (dashed lines). The DW-to-MI/SF boundary (filled circles) has been 
obtained by analyzing the density-wave order parameter, while the MI-SF boundary (filled squares) 
through the charge gap. The two vertical dotted lines denoted two cuts in the phase diagram which will 
be analyzed in details below. The error bars in the points characterizing the DW-to-MI/SF boundary take 
into account the discretization of the g 2 values that we adopted in our numerical simulations (see the 
text). 


stabilized. Their boundaries have been obtained by means of a finite-size scaling of the 
numerical data, for systems up to L — 300 sites. In our simulations we imposed a cutoff 
photon number in each cavity, such that n[’ ho1 < 3. We also truncated the effective Hilbert 
space dimension to a value m = 80 in all the simulations, except for those shown in Fig. 6 
for the neutral energy gap (see the discussion in Sec. 4.3). We checked that, by increasing 
m and the local fock-space truncation over the photon number, the results concerning the 
charge gap and the DW order parameter do not change on the scales shown here. 

For small photon hopping (t/Q 5) 0.2), by increasing the nonlocal nonlinearity g 2 the 
system exhibits a direct transition from the MI to the DW phase. On the other hand, for 
t/£l <; 0.2, the MI-to-DW transition is mediated by an extended region appearing at inter¬ 
mediate g 2 values, where the system stabilizes into a gapless SF. In the following we are 
going to elucidate our finite-size scaling procedure and how we were able to distinguish 
between the different phases. 

4.1 Boundary between Ml and SF phases 

In the limit of small g 2 and t values, the dominant presence of the on-site interactions 
stabilize the system into a MI phase with exactly one polariton per cavity {h = 1), and where 
the charge energy gap has a finite value. As long as the hopping strength t is progressively 
increased (and for fixed g\, g 2 ), the system eventually enters a SF phase, with a vanishing 
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Figure 3 Analysis of the MI-SF Boundary. Left panel: system-size dependence of the charge gap A E c 
per site, for fixed t/Cl = 0.25 and two values of g 2 in the Ml (g 2 /Cl = 1.35) and in the SF (g 2 /Cl = 1.5) 
phase. Symbols denote the DMRG results. Solid and dashed lines are linear and quadratic fitting 
curves, respectively. The difference between the extrapolated values in the two fits AE” = A 

is negligible. Right panel: determination of the critical g* 2 value for the quantum phase transition. The 
triangles denote the charge gap per site A£“ at the thermodynamic limit, as extrapolated in the left 
panel. The dashed line is a best linear fit of the data vs. g 2 . The critical point is obtained when A E c 
vanishes. For t/Cl = 0.25, we get g\/Cl « 1.379. 


gap. The filled squares of Fig. 2 denoting the MI/SF boundaries have been obtained by 
means of a finite-size scaling of the charge gap. We performed simulations up to L = 100 
sites and analyzed whether the gap closes or remains finite in the thermodynamic limit 
L — > oo. 

In Fig. 3, left panel, we highlight the size-dependence of A E c as a function of \/L for 
two points in the phase space close to the MI/SF transition (see points along the dotted 
line in Fig. 2). We expect to see a quadratic dependence A E c ~ L 1 (dashed line) at large 
L [52,53], however a linear extrapolation (solid line) is already a good approximation to 
the scaling, and we can use it to determine A E c in the thermodynamics limit. Indeed, we 
observe that the difference between quadratic and linear extrapolation is tiny (< 10 ') and 
does not produce any distinguishable modification on the scale of Fig. 2. In the specific 
case of Fig. 3, we fixed t/Sl = 0.25 and chose two different values of g 2 /Q corresponding 
to configurations in the gapped MI (g 2 /Cl = 1.35, triangles) and in the gapless SF phase 
{g 2 /£l = 1-5, squares). The MI is signaled by an extrapolated finite value of [ini/_ >co A E c > 
0, while in the SF this is zero. 

In order to locate the critical g 2 for a given value of t (filled squares in Fig. 2) we perform a 
linear extrapolation of the charge gaps in the vicinity of the critical value of g 2 . An example 
of such procedure is shown in the right panel of Fig. 3, where we plot A E c as a function 
of g 2 , when this is close to the phase transition (in the specific, here we set t/Cl = 0.25). 
After a linear extrapolation, we get a critical g 2 value corresponding to g 2 /Q ~ 1.379. An 
analogous procedure is repeated for all the filled squares shown in Fig. 2, thus identifying 
the MI/SF boundary. 

4.2 Boundary of the DW phase 

The DW phase is characterized by a finite order parameter <9 dw- Let us therefore look at 
the two-point staggered correlator in Eq. (4). Since in DMRG simulations we are employ¬ 
ing open boundary conditions, to minimize the border effects we analyze the correlations 
in such a way that the two points are taken symmetrically with respect to the center of the 
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Figure 4 Determination of the DW boundaries. The two-point correlation function CdwM for the 
polariton number and its asymptotic value near the MI-DW quantum phase transition. Here we fix 
t/Q. = 0.05 and vary g2/0.. Left panel: behavior at fixed system size L = 300, as a function of the 
distance r and for different values of g 2 /Ci = 1-3, 1.35, 1.4, and 1.6. To minimize boundary effects, we 
chose the two points (i, i + r) symmetrically with respect to the center of the array. Right panel: finite-size 
scaling close to the transition. Empty circles, filled circles, and triangles respectively are for g 2 /Cl = 1.4 
(DW), 1.35 (near the critical point), and 1.3 (Ml phase). In the Ml phase, Cow vanishes exponentially with 
L, according to the fits: C® 2 ^ 135 « 0.241 x e ~ om2L and C® 2 ^ 13 « 0.253 x e -° 065L , | n the DW phase, Cdw 
converges to a finite value. The inset shows such obtained asymptotic value <5 dw, as a function of g 2 - 


system |21 . The left panel of Fig. 4 shows how differently such polariton correlations behave 
when the system goes from the MI to DW phase, for a fixed system size. 

To be more accurate, in the right panel we performed a finite-size scaling and showed that 
the staggered correlation Cdw(/) approaches the zero value exponentially with L, in the MI 
phase (a similar behavior occurs in the SF region). On the other hand, in the DW such 
correlator asymptotically converges to a finite value. In the specific, here we fix t/Q. = 0.05 
and show that for gi/Q = 1-3, 1.35 the DW order is exponentially suppressed with L, while 
for gi/Q = 1.4 it remains finite. The <9dw order parameter reached for L —> oo is displayed 
in the inset as a function of g2- 

In order to determine the DW boundary in the phase diagram of Fig. 2, we adopted the 
following protocol. For a fixed value of t/Ll, we start increasing g2 from zero up to a finite 
value, with a fixed increment 6 g 2 = 0.05 £2, and to compute the DW order parameter for all 
such values of g2 . The boundary of DW phase in the go —t plane (filled circles in Fig. 2), for 
any fixed f, is located by the g*(f) that gives the first non-vanishing order parameter t? DW . 

Here we stress that, because of the arrangement of our ID array [see Fig. 1(a)] and of 
the asymmetric coupling between the atom and its right/left cavity, the antiferromagnetic 
symmetry of the system is spontaneously broken. In particular, the state |4) of the L-th atom 
will be never occupied, since the transition |2) <-» |4) does not couple to any cavity mode 
[see Eq. (1)]. As a consequence, in our simulations we do not need any symmetry-breaking 
potential. We can observe that the expectation value for the onsite number of polaritons 
explicitly exhibits a staggered behavior, in that the occupation of the (2 n - l)-th site is 
always higher than that of the (2«)-th site (for any integer value of n). Finally we notice 
that such staggering persists at finite size, also for the set of parameters corresponding to 
the MI phase, although it is extremely tiny and decreases with L. This effect eventually 
disappears in the thermodynamic limit. 

[2I The two points of (Snf'srff), with \i-j\ = r, have been chosen such that; = (L-r+i)/ 2 , j = (L+r+i )/2 
for odd r, and ; = (i, - r )/2,j = (L + r)/2 for even r (e.g. for l = too sites, r = l corresponds to 
i = 5o, j = 5i; r = 2 corresponds to / = 49, j = 51; r = 3 to ; = 49 j = 52, and so on) 
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Figure 5 DW phase order parameter as a function of detuning. For a positive detuning, the system will 
never present as a DW phase, meaning that Cdw(L/2) = 0. Taking a negative value of the detuning, we 
observe that Cdw(L/2 ) increases with |A|. The parameters in this figure are t/£l = 0.05, gi/fi = 0.8, 
g 2 /Cl = 1.6 and the system size is L = 100. 


The extension of the DW phase depends on the cavity detuning A. In particular, the ro¬ 
bustness of the order parameter increases with increasing the modulus of the detuning (see 
Fig. 5). Quite remarkably, we note that a positive A will never stabilize an antiferromagnetic 
DW ordering. 

4.3 Neutral gap 

The analysis leading to the phase diagram in Fig. 2 has been corroborated by a study of 
the neutral gap, which vanishes both in proximity of the phase transitions and in the entire 
superfluid region. Differently for the charge gap, it is able to detect the presence of excita¬ 
tions other than particle-hole, and thus locates the boundaries of insulating regions (as the 
DW) beyond the MI. 

The data displayed in Fig. 6 show the behavior of A E n as a function of g 2 , for a fixed 
value of t/Q. In particular we analyzed a vertical cut in the phase diagram of Fig. 2 (see 
the rightmost vertical dotted line in that figure), where the system can be in three different 
phases according to the value of g 2 . With increasing g 2 , it goes from the MI phase (nonzero 
A E„, for t/Q < 1.45) to the SF phase (zero A E n , for 1.45 < t/Q < 1.8), and then to the DW 
phase (nonzero A E n , for t/Q > 1.8). 
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While we cannot see a clear signature of a finite gap for g 2 = 1.8Q, the scaling with the 
size displayed in the left panel of Fig. 6 seems to suggest the scenario depicted above. It is 
however important to stress that the DMRG simulations needed to compute the neutral gap 
have to target the two lowest-lying eigenstates in a single run. Thus they generally require a 
larger dimension m of the effective Hilbert space, as compared to all the other ground-state 
calculations discussed before. The analysis of the neutral gap requires a careful conver¬ 
gence test of the results with m, which we provide in the right panel of Fig. 6. We observe 
that the non monotonic features that are visible in the region 1.45 < t/£l <1.8 have to be 
probably ascribed to the inaccuracy of the method at small m values. This signals the pres¬ 
ence of the gapless SF phase there, in agreement with the results provided by the charge 
gap (MI/SF boundary) and for the DW order parameter (SF/DW boundary). 

5 Summary 

Using the density-matrix renormalization group with open boundary conditions, we studied 
the equilibrium phase diagram of a system of coupled QED cavities in one dimension. We 
provided results beyond the standard model of couplings through photon hopping, and also 
considered nearest-neighbor cross-Kerr nonlinearities. Our analysis is based on a finite-size 
scaling of the ground-state charge and neutral gaps, as well as of the density-wave order 
parameter, for systems up to 300 sites. We showed that, beyond the conventional Mott 
insulator and superfluid phases, the presence of a nearest-neighbor nonlinear coupling can 
also stabilize a density-wave ordering of polaritons. 


Competing interests 

The authors declare that they have no competing interests. 

Authors contributions 

All the authors participated in the design of the research, analysis of the results, and writing of the paper. The DMRG code 
used to run all the simulations of this research has been developed and written by D. Rossini and coworkers [44] (see also 
www.dmrg. it). The DMRG simulations were performed by J. Jin. 

Acknowledgements 

We would like to acknowledge our previous collaboration with M. Hartmann and M. Leib which was inspiring for the present 
work. This work was supported by Italian MIUR via FIRB Project RBFR12NLNA and PRIN Project 2010LLKJBX, by EU 
through IP-SIQS, and by National Natural Science Foundation of China under Grant No. 11175033 and No. 11305021. 

Author details 

1 NEST, Scuola Normale Superiore and Istituto di Nanoscienze - CNR, Piazza dei Cavalieri, 7, 56126 Pisa, Italy. 2 Center 
for Quantum Technologies, National University of Singapore, 117543 Singapore. 

References 

1. Bloch, I., Dalibard, J., Zwerger, W.: Many-body physics with ultracold gases. Rev. Mod. Phys. 80, 885 (2008) 

2. Fisher, M.P.A., Weichman, P.B., Grinstein, G., Fisher, D.S.: Boson localization and the superfluid-insulator transition. 
Phys. Rev. B 40, 546 (1989) 

3. Jaksch, D., Bruder, C., Cirac, J.I., Gardiner, C.W., Zoller, P.: Cold bosonic atoms in optical lattices. Phys. Rev. Lett. 81, 
3108 (1998) 

4. Greiner, M., Mandel, O., Esslinger, T., Hansch, T.W., Bloch, I.: Quantum phase transition from a superfuid to a Mott 
insulator in a gas of ultracold atoms. Nature 415, 39 (2002) 

5. Lewenstein, M., Sanpera, A., Ahufinger, V., Damski. B., Sen(De), A., Sen, U.: Ultracold atomic gases in optical lattices: 
mimicking condensed matter physics and beyond. Adv. Phys. 56, 243 (2007) 

6. Lahaye, T., Menotti, C., Santos, L., Lewenstein, M., Pfau, T.: The physics of dipolar bosonic quantum gases. Rep. Prog. 
Phys. 72, 126401 (2009) 

7. Schmidt, H., Imamoglu, A.: Giant Kerr nonlinearities obtained by electromagnetically induced transparency. Opt. Lett. 
21,1936 (1996) 

8. Imamoglu, A., Schmidt, H., Woods, G., Deutsch, M.: Strongly Interacting Photons in a Nonlinear Cavity. Phys. Rev. Lett. 
79, 1467(1997) 

9. Hartmann, M.J., Brandao, F.G.S.L., Plenio, M.B.: Strongly interacting polaritons in coupled arrays of cavities. Nature 
Phys. 2, 849 (2006) 

10. Greentree, A.D., Tahan, C., Cole, J.H.. Hollenberg, L.C.L.: Quantum phase transitions of light. Nature Phys. 2, 856 
(2006) 


Jin et al. 


Page 11 of 12 


11. Angelakis, D.G., Santos, M.F., Bose, S.: Photon-blockade-induced Mott transitions and XY spin models in coupled 
cavity arrays. Phys. Rev. A 76, 031805(R) (2007) 

12. Hartmann, M.J., Brandao, F.G.S.L., Plenio, M.B.: Quantum many-body phenomena in coupled cavity arrays. Laser & 
Photon. Rev. 2, 527 (2008) 

13. Tomadin, A., Fazio, R.: Many-body phenomena in QED-cavity arrays. J. Opt. Soc. Am. 27, A130 (2010) 

14. Houck, A.A., Tureci, H.E., Koch, J.: On-chip quantum simulation with superconducting circuits. Nature Phys. 8, 292 
( 2012 ) 

15. Schmidt, S., Koch, J.: Circuit QED lattices: Towards quantum simulation with superconducting circuits. Ann. Phys. 525, 
395 (2013) 

16. Carusotto, I., Gerace, D., De Liberato, S., Ciuti, C., Imamoglu, A.: Fermionized photons in an array of driven dissipative 
nonlinear cavities. Phys. Rev. Lett. 103, 033601 (2009) 

17. Tomadin, A., Giovannetti, V., Fazio, R., Gerace, D., Carusotto, I., Tureci, H. E., Imamoglu, A.: Signatures of the 
superfluid-insulator phase transition in laser-driven dissipative nonlinear cavity arrays. Phys. Rev. A 81 ,061801 (R) 
( 2010 ) 

18. Hartmann, M. J.: Polariton crystallization in driven arrays of lossy nonlinear resonators. Phys. Rev. Lett. 104, 113601 
( 2010 ) 

19. Nunnenkamp, A., Koch, J., Girvin, S. M.: Synthetic gauge fields and homodyne transmission in Jaynes-Cummings 
lattices. New J. Phys. 13, 095008 (2011) 

20. Nissen, F., Schmidt, S., Biondi, M., Blatter, G., Tureci, H. E., Keeling, J.: Nonequilibrium dynamics of coupled 
qubit-cavity arrays. Phys. Rev. Lett. 108, 233603 (2012) 

21. Grujic, T., Clark, S. R., Angelakis, D. G., Jaksch, D.: Non-equilibrium many-body effects in driven nonlinear resonator 
arrays. New J. Phys. 14, 103025 (2012) 

22. Grujic, T., Clark, S. R., Jaksch, D., Angelakis, D. G.: Repulsively induced photon superbunching in driven resonator 
arrays. Phys. Rev. A 87, 053846 (2013) 

23. Jin, J., Rossini, D., Fazio, R., Leib, M., Hartmann, M.J.: Photon Solid Phases in Driven Arrays of Nonlinearly Coupled 
Cavities. Phys. Rev. Lett. 110, 163605 (2013) 

24. Jin, J., Rossini, D., Leib, M., Hartmann, M.J., Fazio, R.: Steady-state phase diagram of a driven QED-cavity array with 
cross-Kerr nonlinearities. Phys. Rev. A 90, 023827 (2014) 

25. Underwood, D.L., Shanks, W.E., Koch, J., Houck, A.A.: Low-disorder microwave cavity lattices for quantum simulation 
with photons. Phys. Rev. A 86, 023837 (2012) 

26. Abbarchi, M., Amo, A., Sala, V.G., Solnyshkov, D.D., Flayac, H., Ferrier, L., Sagnes, I., Galopin, E., Lemaitre, A., 
Malpuech, G., Bloch, J.: Macroscopic quantum self-trapping and Josephson oscillations of exciton polaritons. Nature 
Phys. 9, 275 (2013) 

27. Toyoda, K., Matsuno, Y., Noguchi, A., Haze, S., Urabe, S.: Experimental Realization of a Quantum Phase Transition of 
Polaritonic Excitations. Phys. Rev. Lett. Ill, 160501 (2013) 

28. Lucero, E., Barends, R., Chen, Y., Kelly, J., Mariantoni, M., Megrant, A., O’Malley, P, Sank, D., Vainsencher, A., 
Wenner, J., White, T., Yin, Y., Cleland, A.N., Martinis, J.M.: Computing prime factors with a Josephson phase qubit 
quantum processor. Nature Phys. 8, 719 (2012) 

29. Steffen, L., Salathe, Y., Oppliger, M., Kurpiers, R, Baur, M., Lang, C., Eichler, C., Puebla-Hellmann, G., Fedorov, A., 
Wallraff, A.: Deterministic quantum teleportation with feed-forward in a solid state system. Nature 500, 319 (2013) 

30. Chen, Y., Roushan, P, Sank, D., Neill, C., Lucero, E., Mariantoni, M., Barends, R., Chiaro, B., Kelly, J., Megrant, A., 
Mutus, J.Y., O’Malley, P.J.J., Vainsencher, A., Wenner, J., White, T.C., Yin, Y., Cleland, A.N., Martinis, J.M.: Emulating 
weak localization using a solid-state quantum circuit. Nature Commun. 5, 5184 (2014) 

31. Zueco, D., Mazo, J.J., Solano, E., Garcia Ripoll, J.J.: Microwave photonics with Josephson junction arrays: Negative 
refraction index and entanglement through disorder. Phys. Rev. B 86, 024503 (2012) 

32. Peropadre, B., Zueco, D., Wulschner, F., Deppe, F., Marx, A., Gross, R., Garcia Ripoll, J.J.: Tunable coupling 
engineering between superconducting resonators: From sidebands to effective gauge fields. Phys. Rev. B 87, 134504 
(2013) 

33. Sowinski, T., Dutta, O., Hauke, P, Tagliacozzo, L., Lewenstein, M.: Dipolar Molecules in Optical Lattices. Phys. Rev. 
Lett. 108, 115301 (2012) 

34. Dalla Torre, E.G., Berg, E., Altman, E.: Hidden Order in 1D Bose Insulators. Phys. Rev. Lett. 97, 260401 (2006) 

35. Berg, E., Dalla Torre, E.G., Giamarchi, T., Altman, E.: Rise and fall of hidden string order of lattice bosons . Phys. Rev. 
B 77, 245119 (2008) 

36. Rossini, D., Fazio, R.: Phase diagram of the extended Bose-Hubbard model. New J. Phys. 14, 065012 (2012) 

37. Deng, X., Citro, R., Orignac, E., Minguzzi, A., Santos, L.: Polar bosons in one-dimensional disordered optical lattices. 
Phys. Rev. B 87, 195101 (2013) 

38. Wikberg, E-, Larson, J., Bergholtz, E.J., Karlhede, A.: Fractional domain walls from on-site softening in dipolar bosons. 
Phys. Rev. A 85, 033607 (2012) 

39. Batrouni, G.G., Scalettar, R.T., Rousseau, V.G., Gremaud, B.: Competing Supersolid and Haldane Insulator Phases in 
the Extended One-Dimensional Bosonic Hubbard Model. Phys. Rev. Lett. 110, 265303 (2013) 

40. Batrouni, G.G., Rousseau, V.G., Scalettar, R.T., Gremaud, B.: Competing phases, phase separation, and coexistence 
in the extended one-dimensional bosonic Hubbard model. Phys. Rev. B 90, 205123 (2014) 

41. Schiro, M., Bordyuh, M., Oztop, B., Tureci, H.E.: Phase Transition of Light in Cavity QED Lattices. Phys. Rev. Lett. 

109, 053601 (2012) 

42. Kumar, B., Jalal, S.: Quantum Ising dynamics and Majorana-like edge modes in the Rabi lattice model. Phys. Rev. A 
88, 011802(R) (2013) 

43. Schollwock, U.: The density-matrix renormalization group. Rev. Mod. Phys. 77, 259 (2005) 

44. De Chiara, G., Rizzi, M., Rossini, D., Montangero, S: Density Matrix Renormalization Group for Dummies. J. Comput. 
Theor. Nanosci. 5, 1 (2007) 

45. Rossini, D., Fazio, R.: Mott-Insulating and Glassy Phases of Polaritons in 1D Arrays of Coupled Cavities. Phys. Rev. 
Lett. 99, 186401 (2007) 

46. Rossini, D., Santoro, G.E., Fazio, R.: Photon and polariton fluctuations in arrays of QED-cavities. Europhys. Lett. 83, 


Jin et al. 


Page 12 of 12 


47011 (2008) 

47. D’Souza, A.G., Sanders, B.C., Feder, D.L.: Fermionized photons in the ground state of one-dimensional coupled 
cavities. Phys. Rev. A 88, 063801 (2013) 

48. Flu, Y., Ge, G.Q., Chen, S., Yang, X.F., Chen, Y.L.: Cross-Kerr-effect induced by coupled Josephson qubits in circuit 
quantum electrodynamics. Phys. Rev. A 84, 012329 (2011) 

49. Rebic, S., Twamley, J., Milburn, G. J.: Giant Kerr Nonlinearities in Circuit Quantum Electrodynamics. Phys. Rev. Lett. 
103, 150503 (2009) 

50. Verstraete, F., Garcia-Ripoll, J.J., Cirac, J.I.: Matrix Product Density Operators: Simulation of Finite-Temperature and 
Dissipative Systems. Phys. Rev. Lett. 93, 207204 (2004) 

51. Zwolak, M., Vidal, G.: Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent 
Superoperator Renormalization Algorithm. Phys. Rev. Lett. 93, 207205 (2004) 

52. Kiihner, T.D., Monien. H.: Phases of the one-dimensional Bose-Flubbard model. Phys. Rev. B 58, R14741 (1998) 

53. Kuhner, T.D., White, S.R., Monien, H.: One-dimensional Bose-Hubbard model with nearest-neighbor interaction. Phys. 
Rev. B 61, 12474 (2000) 


